home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpkdrv.zip / ST.FOR < prev    next >
Text File  |  1984-01-12  |  16KB  |  546 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL STRTS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE STRTS(LUNIT)
  14.       INTEGER LUNIT
  15. C     LUNIT IS THE OUTPUT UNIT NUMBER
  16. C
  17. C     TESTS
  18. C        STRCO, STRSL
  19. C
  20. C     LINPACK. THIS VERSION DATED 08/14/78 .
  21. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  22. C
  23. C     SUBROUTINES AND FUNCTIONS
  24. C
  25. C     LINPACK STRCO,STRSL
  26. C     EXTERNAL STRXX,SMACH
  27. C     BLAS SAXPY,SDOT,SSCAL,SASUM
  28. C     FORTRAN ABS,AMAX1,FLOAT,MAX0
  29. C
  30. C     INTERNAL VARIABLES
  31. C
  32.       REAL A(15,15),B(15),BT(15),X(15),XEXACT(15),XT(15),Z(15)
  33.       REAL AINV(15,15),DET(2)
  34.       REAL SDOT,STUFF,T
  35.       REAL ANORM,AINORM,RCOND,COND,COND1,SMACH,EPS
  36.       REAL ENORM,ETNORM,RNORM,RTNORM,XNORM,XTNORM,EN,SASUM
  37.       REAL FNORM,ONEPX,Q(7),QS(7)
  38.       INTEGER I,INFO,J,JOB,KASE,KOUNT,KSING,LDA,ML,MU,N,NPRINT
  39.       INTEGER KSUSP(7),IQ(7)
  40. C
  41.       LDA = 15
  42. C
  43. C     WRITE MATRIX AND SOLUTIONS IF  N .LE. NPRINT
  44. C
  45.       NPRINT = 3
  46. C
  47.       WRITE (LUNIT,380)
  48.       WRITE (LUNIT,730)
  49. C
  50.       DO 10 I = 1, 7
  51.          KSUSP(I) = 0
  52.    10 CONTINUE
  53.       KSING = 0
  54. C
  55. C     SET EPS TO ROUNDING UNIT
  56. C
  57.       EPS = SMACH(1)
  58.       WRITE (LUNIT,390) EPS
  59.       WRITE (LUNIT,370)
  60. C
  61. C     START MAIN LOOP
  62. C
  63.       KASE = 1
  64.    20 CONTINUE
  65. C
  66. C        GENERATE TEST MATRIX
  67. C
  68.          CALL STRXX(A,LDA,N,KASE,LUNIT)
  69. C
  70. C        N = 0 SIGNALS NO MORE TEST MATRICES
  71. C
  72. C     ...EXIT
  73.          IF (N .LE. 0) GO TO 360
  74.          ANORM = 0.0E0
  75.          DO 30 J = 1, N
  76.             ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
  77.    30    CONTINUE
  78.          WRITE (LUNIT,570) ANORM
  79. C
  80.          IF (N .GT. NPRINT) GO TO 50
  81.             WRITE (LUNIT,370)
  82.             DO 40 I = 1, N
  83.                WRITE (LUNIT,600) (A(I,J), J = 1, N)
  84.    40       CONTINUE
  85.             WRITE (LUNIT,370)
  86.    50    CONTINUE
  87. C
  88. C        GENERATE EXACT SOLUTION
  89. C
  90.          XEXACT(1) = 1.0E0
  91.          IF (N .GE. 2) XEXACT(2) = 0.0E0
  92.          IF (N .LE. 2) GO TO 70
  93.             DO 60 I = 3, N
  94.                XEXACT(I) = -XEXACT(I-2)
  95.    60       CONTINUE
  96.    70    CONTINUE
  97. C
  98. C        GENERATE R.H.S.
  99. C
  100.          DO 90 I = 1, N
  101.             B(I) = 0.0E0
  102.             BT(I) = 0.0E0
  103.             DO 80 J = 1, N
  104.                B(I) = B(I) + A(I,J)*XEXACT(J)
  105.                BT(I) = BT(I) + A(J,I)*XEXACT(J)
  106.    80       CONTINUE
  107.             X(I) = B(I)
  108.             XT(I) = BT(I)
  109.    90    CONTINUE
  110. C
  111. C        UPPER OR LOWER TRIANGULAR
  112. C
  113.          ML = 0
  114.          MU = 0
  115.          DO 120 J = 1, N
  116.             DO 110 I = 1, N
  117.                IF (A(I,J) .EQ. 0.0E0) GO TO 100
  118.                   IF (I .LT. J) MU = MAX0(MU,J-I)
  119.                   IF (I .GT. J) ML = MAX0(ML,I-J)
  120.   100          CONTINUE
  121.   110       CONTINUE
  122.   120    CONTINUE
  123.          WRITE (LUNIT,670) ML,MU
  124.          IF (ML .NE. 0 .AND. MU .NE. 0) GO TO 350
  125.             IF (MU .EQ. 0) JOB = 0
  126.             IF (ML .EQ. 0) JOB = 1
  127.             IF (JOB .EQ. 0) WRITE (LUNIT,710)
  128.             IF (JOB .EQ. 1) WRITE (LUNIT,720)
  129.             STUFF = 4095.0E0
  130.             DO 140 J = 1, N
  131.                DO 130 I = 1, N
  132.                   IF (I .LT. J .AND. JOB .EQ. 0) A(I,J) = STUFF
  133.                   IF (I .GT. J .AND. JOB .EQ. 1) A(I,J) = STUFF
  134.   130          CONTINUE
  135.   140       CONTINUE
  136. C
  137. C           ESTIMATE CONDITION
  138. C
  139.             CALL STRCO(A,LDA,N,RCOND,Z,JOB)
  140. C
  141. C           OUTPUT NULL VECTOR IF N .LE. NPRINT
  142. C
  143.             IF (N .GT. NPRINT) GO TO 160
  144.                WRITE (LUNIT,620)
  145.                DO 150 I = 1, N
  146.                   WRITE (LUNIT,630) Z(I)
  147.   150          CONTINUE
  148.                WRITE (LUNIT,370)
  149.   160       CONTINUE
  150. C
  151. C
  152. C           TEST FOR SINGULARITY
  153. C
  154.             IF (RCOND .GT. 0.0E0) GO TO 170
  155.                WRITE (LUNIT,610) RCOND
  156.                WRITE (LUNIT,400)
  157.                KSING = KSING + 1
  158.             GO TO 340
  159.   170       CONTINUE
  160.                COND = 1.0E0/RCOND
  161.                WRITE (LUNIT,420) COND
  162.                ONEPX = 1.0E0 + RCOND
  163.                IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,410)
  164. C
  165. C              COMPUTE INVERSE, DETERMINANT AND COND1 = TRUE CONDITION
  166. C
  167.                DO 190 J = 1, N
  168.                   DO 180 I = 1, N
  169.                      AINV(I,J) = A(I,J)
  170.   180             CONTINUE
  171.   190          CONTINUE
  172.                CALL STRDI(AINV,LDA,N,DET,110+JOB,INFO)
  173.                AINORM = 0.0E0
  174.                DO 200 J = 1, N
  175.                   IF (JOB .EQ. 0)
  176.      *               AINORM = AMAX1(AINORM,SASUM(N-J+1,AINV(J,J),1))
  177.                   IF (JOB .EQ. 1)
  178.      *               AINORM = AMAX1(AINORM,SASUM(J,AINV(1,J),1))
  179.   200          CONTINUE
  180.                COND1 = ANORM*AINORM
  181.                WRITE (LUNIT,430) COND1
  182.                WRITE (LUNIT,650) DET(1)
  183.                WRITE (LUNIT,660) DET(2)
  184. C
  185. C              SOLVE  A*X = B  AND  TRANS(A)*XT = BT
  186. C
  187.                CALL STRSL(A,LDA,N,X,JOB,INFO)
  188.                CALL STRSL(A,LDA,N,XT,JOB+10,INFO)
  189. C
  190.                IF (N .GT. NPRINT) GO TO 230
  191.                   WRITE (LUNIT,440)
  192.                   DO 210 I = 1, N
  193.                      WRITE (LUNIT,640) X(I)
  194.   210             CONTINUE
  195.                   WRITE (LUNIT,450)
  196.                   DO 220 I = 1, N
  197.                      WRITE (LUNIT,640) XT(I)
  198.   220             CONTINUE
  199.                   WRITE (LUNIT,370)
  200.   230          CONTINUE
  201. C
  202. C              RESTORE ZEROS IN OTHER TRIANGLE
  203. C
  204.                DO 260 J = 1, N
  205.                   DO 250 I = 1, N
  206.                      IF (A(I,J) .NE. STUFF) GO TO 240
  207.                         A(I,J) = 0.0E0
  208.                         AINV(I,J) = 0.0E0
  209.   240                CONTINUE
  210.   250             CONTINUE
  211.   260          CONTINUE
  212. C
  213. C              COMPUTE ERRORS AND RESIDUALS
  214. C                 E  =  X - XEXACT
  215. C                 ET =  XT - XEXACT
  216. C                 R  =  B - A*X
  217. C                 RT =  BT - A*XT
  218. C                 AI = A*INV(A) - I
  219. C
  220.                XNORM = SASUM(N,X,1)
  221.                XTNORM = SASUM(N,XT,1)
  222.                ENORM = 0.0E0
  223.                ETNORM = 0.0E0
  224.                FNORM = 0.0E0
  225.                DO 270 J = 1, N
  226.                   ENORM = ENORM + ABS(X(J)-XEXACT(J))
  227.                   ETNORM = ETNORM + ABS(XT(J)-XEXACT(J))
  228.                   T = -X(J)
  229.                   CALL SAXPY(N,T,A(1,J),1,B,1)
  230.                   BT(J) = BT(J) - SDOT(N,A(1,J),1,XT,1)
  231.   270          CONTINUE
  232.                RNORM = SASUM(N,B,1)
  233.                RTNORM = SASUM(N,BT,1)
  234. C
  235. C
  236. C              A*INV(A) - I
  237. C
  238.                AINORM = 0.0E0
  239.                DO 300 J = 1, N
  240.                   DO 280 I = 1, N
  241.                      B(I) = 0.0E0
  242.   280             CONTINUE
  243.                   DO 290 K = 1, N
  244.                      T = AINV(K,J)
  245.                      CALL SAXPY(N,T,A(1,K),1,B,1)
  246.   290             CONTINUE
  247.                   B(J) = B(J) - 1.0E0
  248.                   AINORM = AMAX1(AINORM,SASUM(N,B,1))
  249.   300          CONTINUE
  250. C
  251.                WRITE (LUNIT,460) ENORM,ETNORM
  252.                WRITE (LUNIT,470) RNORM,RTNORM
  253.                WRITE (LUNIT,580) AINORM
  254. C
  255. C              COMPUTE TEST RATIOS
  256. C
  257.                Q(1) = COND/COND1
  258.                Q(2) = COND1/COND
  259.                Q(3) = ENORM/(EPS*COND*XNORM)
  260.                Q(4) = ETNORM/(EPS*COND*XTNORM)
  261.                Q(5) = RNORM/(EPS*ANORM*XNORM)
  262.                Q(6) = RTNORM/(EPS*ANORM*XTNORM)
  263.                Q(7) = AINORM/(EPS*COND)
  264.                WRITE (LUNIT,370)
  265.                WRITE (LUNIT,480)
  266.                WRITE (LUNIT,370)
  267.                WRITE (LUNIT,540)
  268.                WRITE (LUNIT,550)
  269.                WRITE (LUNIT,560)
  270.                WRITE (LUNIT,370)
  271.                WRITE (LUNIT,590) (Q(I), I = 1, 7)
  272.                WRITE (LUNIT,370)
  273. C
  274. C              LOOK FOR SUSPICIOUS RATIOS
  275. C
  276.                QS(1) = 1.0E0 + 4.0E0*EPS
  277.                QS(2) = 10.0E0
  278.                EN = FLOAT(N)
  279.                IF (N .EQ. 1) EN = 2.0E0
  280.                DO 310 I = 3, 7
  281.                   QS(I) = EN
  282.   310          CONTINUE
  283.                KOUNT